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ABSTRACT 

We model fluctuations in the Cosmic Infrared Background (CIB) arising from known galaxy popu- 
lations using 233 measured UV, optical and NIR luminosity functions (LF) from a variety of surveys 
spanning a wide range of redshifts. We compare best-fit Schechter parameters across the literature 
and find clear indication of evolution with redshift. Providing fitting formulae for the multi-band evo- 
lution of the LFs out to z^5, we calculate the total emission redshifted into the near-IR bands in the 
observer frame and recover the observed optical and near-IR galaxy counts to a good accuracy. Our 
empirical approach, in conjunction with a halo model describing the clustering of galaxies, allows us 
to compute the fluctuations of the unresolved CIB and compare the models to current measurements. 
We find that fluctuations from known galaxy populations are unable to account for the large scale 
CIB clustering signal seen by Spitzer/IRAC and AKARI/IRC and continue to diverge out to larger 
angular scales. This holds true even if the LFs are extrapolated out to faint magnitudes with a steep 
faint-end slope all the way to z—8. We also show that removing resolved sources to progressively 
fainter magnitude limits, isolates CIB fluctuations to increasingly higher redshifts. Our empirical 
approach suggests that known galaxy populations are not responsible for the bulk of the fluctuation 
signal seen in the measurements and favors a very faint population of highly clustered sources. 

Subject headings: cosmology: diffuse radiation — large-scale structure of universe — galaxies: evolu- 
tion — luminosity function — infrared radiation 



1. INTRODUCTION 

Cosmic infrared background (CIB) includes contribu- 
tions from emissions over the entire history of the Uni- 
verse, including from objects inaccessible to the cur- 
rent telescopic studies. Several direct measurements of 
the total mean levels of the CIB using the wide-beam 
DIRBE and IRTS instruments claim a significant ex- 
cess mean flux ov er the contribution of known galax- 
ies in the near- IR (iDwek fc ArendtlH998l: iGoriian et all 
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20001: IWright fc Reesel 120001: iCambresv et at. _ 
Matsumoto et all 12005( 1: also see review by IKash linskv 
(|2005l ). The entire excess emission over that from 
know n galaxy populatio ns (~ 30 nWm _2 sr _1 in 1- 
4/im, Kashlinsky (2005)) was orig inally theorized to 
come from primordial Pop III stars ( jSantos et all [20021 : 
iSalvaterra fc Ferraral l2003| ) but this interpretation has 
been challenged on several grounds since the claimed 
levels r equire uncomfortable levels of star formation eff i- 
ciency (|Madau fc Silkll2005l: ISalvaterra fc Ferrarall2006() . 
It is possible that much of the excess flux seen may be 
due to inaccurat e removal of bright zodiacal emis sion in 
the foreground (|Dwek et all [20051 : iMattilal 120011 . Fur- 
thermore, the extragalactic background light (EBL) is a 
fundamental source of opacity for high energy photons 
and the 7-ray attenuation seen in blaz ar spectra favors 
low levels of NIR backgrou nd light fe.g. lAharonian et all 
l2006t iMazTn fc Raudl2007l) . 

An alternative way to study the CIB, much less sen- 
sitive to foreground removal, is to measure background 
anisotropics after subtracting resolved gal axies down to 
faint magnitudes (IKashlinskv et all 1 19961 ). Detections 
of spatial structure in the CIB were initially based on 



datas ets from CO BE/BIRBE (IKashlinskv fc Odeiiwaldl 
2000) , the IRTS ([Matsumoto et all 120001) and 2MASS 
(IKashlinskv et all hooa IQdenwald et all 120031) . More 
recently. IKashlinskv etaL I (120051 . 1200711 using deep 
exposures from Spitzer/IRAC (3.6-8.0/im) found sig- 
nificant fluctuations after subtracting galaxies down 
to ?tias~25. The level of these fluctuations, ~0.1 
nWm _2 sr _1 at arcminute scales, imply an isotropic CIB 
flux as low as ~1 nWm _2 sr _1 from t he remaining un- 
resolv ed sources in the IRAC ba nds (jKashlinsk v et all 
l2007d ). iThompson et all (|2007aD analysis constrained 
CIB at 1.4-1.8/xm usi ng iJST/NICMOS observations 
and IMatsumoto et al.l (|2011D measure fluctuations 
on arcminute scales in the 2.4-4. l^tm range using 
the AK ARI satellite. After this paper was sub- 
mitted, IKashlinsk v et all (|2012l) have measured the 
Spitzer/IRAC out to <1° using more extensive datasets 
from the Spitzer Extended Deep Survey (SEDS), confirm- 
ing earlier results and extending the fluctuation measure- 
ment to much larger angular scales. All the present mea- 
surements of CIB-fluctuations are consistent with an ex- 
tragalactic origin, necessitating an associated unresolved 
component in the CIB. This component likely requires 
only a fraction of the CIB excess, whi ch is below limits 
imposed by 7-ray photon absorption (|Kashlinsky et all 
[2?Wcl IKashlinskv fc Bandl 120071 : lArendt et al.H201f)ll ~ 
There seems to be an emerging consensus that the 
extragalactic clustering signal is real, but the nature 
of the sources producing it is still a subject of de- 
bate. Plausible candidates for the bulk of the CIB are 
evolving stellar populations in galaxies, although ac- 
creting black holes at h igh-z c an also contribut e (e.g. , 
iRicotti fc Ostrikerll200l . Both IKashlinskv et all ( 2005) 
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and iMatsumo to et all (|2011f) argue that that the clus- 
tering is cons is tent w i th "first stars " era o bjects whereas 
iCoorav et ail (|2007t ); iCharv et all (120081 ) have posited 
that the signal originates mostly in the clustering of faint 
galaxies at redshifts z~l-3. Understanding the expected 
levels of fluctuations from known galactic populations is 
possible following the establishment of the standard cos- 
mologic al model for structure formation, the concordance 
ACDM (jKomatsu et all 1201 J) . In order to compute the 
levels of source-s ubtracted CIB fluctu ations remaining in 
the Spitzer data, iSullivan et al.l (|2007l ) used a halo model 
combined with conditional luminosity functions and com- 
pared it to measurements a t 3.6/xm. Their cl a im is that 
the fluctuations detected by Kashlins kv et al.l (|2005l) can 
be explained by ordinary galaxies just beyond the de- 
tection threshold of Spitzer /IRAC, although this claim 
appears to contradict the results of their analysis shown 
i n their Fig. 8 . 

iKashlinskvl (|2005() discusses the importance of the 
shape of the emission history for the resulting fluctua- 
tions demonstrating how brief episodes of light produc- 
tion can lead to enhanced fluctuations. In this paper, we 
construct the entire history of light production produced 
by known galaxy populations using a novel empirical ap- 
proach that relies exclusively on observations. We use a 
compilation of galaxy luminosity functions (LF) in the 
literature to populate the observed lightcone with galax- 
ies down to faint magnitudes. The many galaxy surveys 
conducted in recent years provide a wealth of data in 
multiple bands and cover a wide range of redshifts. Indi- 
vidually, LFs only probe specific rest-frame wavelengths 
for a limited range of redshifts, while together we can 
use them to infer the source distribution composing the 
background light in the 0.1-5.0/im range. Our only the- 
oretical assumptions concern the clustering properties of 
the unresolved sources which are modeled according to 
the well-established concordan c e ACD M model (see Sec- 
tion [5]) . We refer to Uohnstonl ([201 ID for a good review 
on the properties of luminosity functions and how they 
are measured. 

Modeling the underlying populations of the EBL 
has been attempted using various mixtures of theory 
and observations. Backward evolution scenarios take 
the present galaxy popul ations and extrapolate them 
to higher redshift (e.g ., Uimenez fc Kashlinskv! 119991 : 
iFrance schini ct al. 2008), while forward evolution follows 
dark matter merger trees starting from the cosmolog- 
ical initial condit ions, using semi-analytical models o f 
galaxy formation (iGilmore et al.ll20Tol IGuo et alj|2011l ). 
iDommguez et al.l ( 20111 ) u se directly the measu red K- 
band LFs out to z=4 from lCirasuolo et al.l (1—0 1 i. com- 
bined with best-fit SEDs of multi- wavelength galaxy data 
(AEGIS) to empirically derive the overall EBL spec- 
trum. We however, present an alternative empirical ap- 
proach by examining the best-fit Schechter parameters 
(jSchechterl I1976D of 233 LFs covering the UV, optical 
and near-IR out to redshifts z^3-8. We provide empir- 
ical fitting formulae describing the smooth evolution of 
multi-band LFs with redshift, and construct lightconcs 
containing all populations seen in the near-IR bands, se- 
lected at each redshift such that X^ir = 0- + z)\ rest . 

This paper is organized as follows. Section [5] describes 
the data used and in Section [3] we explain the modeling 
in detail. In Section [4j we calculate both galaxy number 



counts and the EBL in the near-IR bands (JHKLM) 
and compare with existing data. In Section [5] we ana- 
lyze the source-subtracted CIB-fluctuations implied by 
our empirical reconstruction and compare with previ- 
ous work. We discuss the implications of our findings 
in Section [6l Throughout this paper we adopt the con- 
cordance ACDM comology with f2 m =0.3, fl\—0.7 and 
Hq=70 km-s -1 -Mpc -1 . All m agnitudes are in th e AB 
system unless stated otherwise (|Oke fc Gunnlll983l) . 

2. MEASUREMENTS OF THE GALAXY 
LUMINOSITY FUNCTION 

The total emission seen in the near-IR bands 
(JHKLM) depends on the contribution of local near-IR 
galaxies as well as redshifted light radiated at shorter 
rest-frame wavelengths. To quantify the present day 
background produced by galaxies, we have utilized mea- 
surements of luminosity functions probing all rest-frame 
wavelengths in the interval 0.1< A <5.0/im anywhere in 
the redshift cone. This results in a compilation of 233 
LFs from a large variety of surveys which we list in Ta- 
ble [TJ Our approach does no t depend on stellar pop- 
ulation synthesis models (e.g.. iBruzual fc Charlotll2003|) 
and we do not need to make an assumption for the IMF. 
Rather, in this method we predict the levels of CIB fluc- 
tuations directly from the available data, assuming only 
i) standard ACDM model of structure formation and ii) , 
the validity of a Schechter-type LF after fitting its pa- 
rameters to the data. All the LFs we use have been 
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Fig. 1. — All 233 luminosity functions used in our analysis in 
Schechter parametrization (see original references in Table [lj. The 
wavelength bins are listed in the panels (lower right) and their effec- 
tive wavelengths are listed in Table [2] along with other properties. 
The LFs shown have a range of redshifts. 



Reference 


Rest-frame band 


Redshift 


Sample 


Selection 


Survey Catalog / Field 


Symbol / Color 6 






z 


N aa l 


m Hm {AB) 


Arnouts ct al. (2005) 


1500A 


0.2-1.2 


1039 


NUV<24.5 


GALEX/VVDS 


green triangles(up) 






1.75-3.4 




F450&F606<27 


HDF 




Wvder et al. (2005) 


NUV,FUV 


0.055 


896,1124 


m uv < 20 


GALEX/2dF 


blue circles 


Oesch et al. (2010) 


1500A 


0.5-2.5 


284-403 


$26 


HST ERS 


yellow circles 


Oesch et al. (20121 


1500 A 


~8 


70 


H<27.5 


CANDLES /HUDF09/ERS 


pink triangles (up) 


Reddv et al. (20081 


1700 A 


1.9-3.4 


~15,000 


7£<25.5 


a 


blue crosses 


Yoshida et al. (2006) 


1500 A 


~4,5 


3808,539 


<26-27 


Subaru Deep Field 


blue squares 


McLure et al. (2009) 


1500 A 


-5,6 


~1500 


z'<26 


SXDS/UKIDSS 


purple squres 


Ouchi et al. (2009) 


1500 A 


7 


22 


<26 


SDF/GOODS-N 




Bouwens et al. (20071 


1600A,1350A 


~4,5,6 


4671,1416,627 


<29 


HUDF/GOODS 


violet triangles (down) 


Bouwens et al. (2011) 


1600A,1750A 


-7,8 


73,59 


<26-29.4 


HUDF09 


orange diamonds 


Gabasch et al. (20041 


u'g> 


0.45-5 


5558 


/<26.8 


FORS Deep Field 


Green triangles(down) 


Baldrv et al. (20051 


°- 1 u 


<0.3 


43223 


u<20.5 


SDSS 


red squares 


Faber et al. (20071 


B 


0.2-1.2 


~34000 


R $ 24 


DEEP2/COMBO-17 


tan squares 


Norbere et al. (2002) 


bj 


<0.2 


110500 


<19.45 


2dFGRS 


purple squares 


Blanton et al. (2003b1 


01 ugriz 


0.1 


147986 


<16.5-18.3 


SDSS 


blue plus 


Montero-Dorta & Prada (2009) 


01 ugriz 


$0.2 


947053 


<17-19 


SDSS 


green crosses 


Lovedav et al. (2012) 


01 ugriz 


0.002-0.5 


8647-12860 


r<19.8 


GAMA 


yellow squares 


Ilbert et al. (2005) 


UBVRI 


0.05-2.0 


11034 


7<24 


VIMOS-VLT Deep Survey 


pink triangles (up) 


Gabasch et al. (2006) 


i'z'r' 


0.45-3.8 


5558 


/<26.8 


FDF 


green circles 


Marchesini et al. (2007) 


BVR 


2.0-3.5 


989 


K s <25 


MUSYC/FIRES /GOODS/EIS 


orange circles 


Marchesini et al. (2012) 


V 


0.4-4.0 


19403 


H<27.8,K"<25.6 


a 


blue triangles(up) 


Hill et al. (2010) 


ugriz 


0.0033-0.1 


2437-3267 


<18-21 


MGC/UKIDSS/SDSS 


purple diamonds 




YJHK 




1589-1798 


<17.5-18 






Dahlen et al. (2005) 


UBR 


0.1-2 


18381 


i?<24.5 


GOODS-HST/CTIO /ESO 


dark green diamonds 




J 


0.1-1 


2768 


K s <23.2 






Jones et al. (2006) 


b 3 r f 


<0.2 


138226 


bjTf <15.6, 16.8 


6dFGS/2MASS 


dark red plus 




JHK 






JHK< 14.7 


/SuperCOSMOS 




Bell et al. (2003) 


ugriz 


< 0.1 


22679 


r<17.5 


SDSS 


orange circles 




K 




6282 


K<15.5 


2MASS 




Kashikawa ct al. (2003) 


BK' 


0.6-3.5 


439 


K'<2i 


Subaru Deep Survey 


red crosses 


Stefanon fc Marchesini (2011) 


JH 


1.5-3.5 


3496 


K s < 22.7-25.5 


MUSYC/FIRES /FIREWORKS 


green squares 


Pozzetti et al. (2003) 


JK S 


0.2-1.3 


489 


i^ s <20 


K20 Survey 


tan plus 


Feulner et al. (2003) 


JK' 


0.1-0.6 


500 


K'< 19.4-20.9 


MUNICS 


yellow crosses 


Eke et al. (2005) 


JK a 


0.01-0.12 


16922,15664 


JK 3 <15.5 


2dFGRS/2MASS 


violet diamonds 


Cole et al. (2001) 


JK S 


0.005-0.2 


7081,5683 


JK S < 15.5 


2dFGRS/2MASS 


blue squares 


Smith et al. (2009) 


K 


0.01-0.3 


40111 


ii"<17.9,r<17.6 


UKIDSS-LAS/SDSS 


red triangle(up) 


Saracco et al. (2006) 


K 8 


0.001-4 


285 


Ks< 24.9 


HDFS/FIRES 


blue triangles (down) 


Kochanek et al. (2001) 


K s 


0.003-0.03 


4192 


K 2Q < 13.35 


2MASS/CfA2/UZC 


magenta circles 


Huane et al. (2003) 


K 


0.001-0.57 


1056 


K<15 


2dF/AAO 


violet diamonds 


Arnouts et al. (2007) 


K 


0.2-2 


21200 


ma.Gmic <21.5 


SWIRE/VVDS 


dark green squares 












/UKIDSS/CFHTLS 




Cirasuolo et al. (2010) 


K 


0.2-4 


~50000 


ii"<23 


UKIDSS/SXDS 


orange plus 


Babbedge et al. (2006) 


i3.6Mm-M4.5/jm 


0.01-0.6 


34281 


<20.2 


SWIRE/INT WFS 


blue crosses 


Dai ct al. (2009) 




0.01-0.6 


4905,5847 


LM <19,I< 20 .4 


IRAC-SS/AGES 


dark red circles 



Note. — The measured LF are shown in Figure [TJ and all Schechter parameters are displayed in Figure [3] The compile database of the Schechter parameters is available upon 
request. 

°Data taken from multiple surveys/fields 

6 The symbols and color of the corresponding data points in Figure [2] and [3] 
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characterized by a Schechter function (lSchechter]|1976l) . 
4>{M)dM = 0.4 In (10)0* (io°- 4 ( M *- M )) a+1 

xexp(-10 - 4(M *- M) )dM, 

(1) 

determined by the normalization, 4>* , characteristic ab- 
solute magnitude, M* and the faint-end slope, a. By 
integrating Equation (p}, the luminosity density can be 
shown to be L — 4>*L*T(a + 2), where L* is the charac- 
teristic luminosity and T(x) is the Gamma function. All 
the Schechter LFs used are shown in Figure [TJ 

The Schechter LF is usually found to fit the data fairly 
accurately but deviations are seen, in particular when 
fitting a wide range of luminosities. At low-z for exam- 
ple, IJones et"al . (2006) find that the sha pe does not fit 
the sha rp do wnturn seen at M* and both [Blanton et al.l 
(|2003bD and IMontero-Dorta fc Pradal (|2009ft find an ex- 
cess of bright galaxies in the blue SDSS bands. There 
are also hints of an upturn in the local LF at faint mag- 
nitudes where the Schechter fit does a poor job (e.g., 
IBlanton et al.l 120051) . We address this faint-end issue in 
Section 13.11 but note that sources at the bright end are 
efficiently removed from the maps in CIB fluctuations 
studies. At longer wavelengths (>5/im), a double power- 
law is found to pr ovide a more adequate fit than the 
Schec hter function (Babbcd ge et al.l l2006; Ma gnelli et al.l 
1201 If) . The mutual consistency of measurements is a pri- 
mary concern when comparing LFs in the literature. In- 
consistencies can be caused by field-to-field variations, 
photometric system, k-corrections, type of LF-estimator, 
survey depth and completeness, redshift binning, sam- 
ple statistics, error estimates, etc. These undoubtedly 
account for differences in shape and amplitude of the 
measured LF (see Figure [3]). We include a discussion of 
common issues in Appendix [Bj but these do not affect 
our results because we let all measurements collectively 
contribute to our derived LF (see Section [3]). 

To directly compare flux measurements at different 
wavelengths, we have adopted the AB magnitude sys- 
tem which conveniently relates the apparent magnitude, 
rriAB, to the specific flux, /„, via 

U = io-°- 4 < m ^- 23 - 9 VJy, (2) 

(1 Jy = 10~ 26 Wm~ 2 Hz _1 ). Where system conver- 
sions are not explicitly given by the authors (Table [lj 
we have made use of the ca lculations available at 
http://mips.as.arizona.edu/^cnaw/sun.html. With all 
magnitudes converted to AB we do not distinguish be- 
tween magnitudes of different filter- and photometric 
variations, e.g. Johnson U and SDSS u, apart from their 
center frequencies. 

3. POPULATING THE LIGHTCONE WITH 
KNOWN GALAXY POPULATIONS 

This section outlines the step-by-step approach lead- 
ing to the quantification of the galaxy distribution seen 
on the sky. Using the data in Table [TJ we populate the 
evolving lightcone by placing the rest-frame galaxy dis- 
tribution at a distance such that the associated emission 
is shifted into the near-IR bands in the observer frame, 
defined by \nir/(1- + z). Initially, we bin the LFs ac- 
cording to their rest-frame wavelength in fiducial bands 




0.1 1.0 10.0 

Wavelength (yu.m) 



Fig. 2. — The local luminosity density according to all available 
LF measurements at z<0.12 in Table [T] with symbols/colors in- 
dicated in the same Table. To avoid overcrowding the region of 
interest we omit error bars. The solid line shows the luminosity 
density in our fiducial bands as implied by our fits in Figure [3] 
The sets of gray lines show the contribution from galaxies of dif- 
ferent metallicities a nd ages from syn thetic galaxy SED spectra 
shown in Fig. 14 of {Kashlinsky 2005). The bottom-gray curves 
show the early type stellar populations, the upper-dark show late 
type populations and middle-light lines show the average of the 
two contributions. 



which we call UV, U, B, V, R, I, z, J, H, K, L and M 
(see Table [2]). For example, measurements in rest- frame 
SDSS g' , Johnson B and 2dF bj are binned together in 
our B-band despite having an offset in center wavelength 
of about 0.03/im. The largest offset occurs in our J-bin 
where the centers of SDSS i and Johnson / is 0.063/im. 
The uncertainty associated with the redshift of the pop- 
ulation usually dominates these offsets so we do not cor- 
rect for them. The centers of our fiducial bands, A e //, 
are taken to be the mean rest-frame wavelength of all 
measurements in the bin (see Table [2]). 

By placing the entire population of each LF at the me- 
dian redshift of the sample, z me d, we examine the evolu- 
tion of the individual Schechter parameters (a, M* ,</)*) in 
our fiducial bands. In the cases where z me d is not explic- 
itly given by the authors, we choose the midpoint of the 
redshift bin of the LF measurement. The distances of the 
galaxies composing the LF is the dominant uncertainty 
in the resulting counts on the sky and we have therefore 
examined the effects of placing the LF at the opposite 
boundaries of the redshift bin (the resulting counts dif- 
fer by less than a factor of two at the two extremes (see 
Section |4|)). Figure |3] shows the Schechter parameters as 
a function of redshift from 0.15-4.5/^m. Across the spec- 
trum, we see clear indication of evolution in M* and 
and in some cases also in the poorly measured a. 

Over time, galaxy populations evolve both in bright- 
ness and abundance. As small systems merge to form 
more massive ones, we expect a net increase in the num- 
ber of bright and massive galaxies with time accom- 
panied by a decrease in fainter ones. This is encoded 
in the evolution of 4>* (the number density of L* sys- 
tems), which we expect to increase with time whereas 



UVUBVR I zJHKLM 




0.1 1 0.1 1 0.1 1 0.1 1 0.1 1 0.1 1 0.1 1 0.1 1 0.1 1 0.1 1 0.1 1 0.1 1 



Redshift 

Fig. 3. — The measured Schechter parameters a,M* ,cf>* from the studies in Table [l] including the luminosity density, C v =cj>* L*T(a + 2), as a function of redshift. The different 
sybmols/colors are listed along with the corresponding references in Table [l] We have omitted error bars for the sake of cl arity. The solid curves show the evolutionary fits according to 
Equations (f3j>— with the best-fit parameters listed in Table[2] We have modified M^ v to follow the fitting functions of Bouwcns ct al. (2011) at z>3.5 to better match the turnover 
seen. We note that our fits are only empirically supported for z<4, beyond which we extrapolate. The dashed curves in the a panels shows the evolution assumed in our default 
model whereas the dark shades area s en compass the range bracketed by our high faint- end (HFE) and low faint-end (LFE) models. These ranges are ultimately constrained by the 
observed galaxy counts (see Section |3. IE ■ The shaded areas in the bottom row (vC v ) is the evolving quantity <f>(z)*vL v (z)*T(ot(z) + 2) corresponding to this allowed range in a(z). 

The dotted curves in (j>* in L and A/-bands are not fits to the data but are instead assumed to have the same form as the K-haad fits. The light gray shaded areas correspond to the 
redshift regions for which the rest-frame emission redshifts into the observed NIR wavelengths of interest, defined to encompass the 1.25-4.5/xm range. We are most concerned with 
the goodness of fit in these regimes. All the data-points assume h=0.7. 
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the faint-end slope, a should consequently flatten. The 
difference of the LF among rest-frame bands reflects the 
tendency of galaxies of different types being preferen- 
tially bright / faint at a given wavelength. The decompo- 
sition of the LF into red/blue galaxies typically shows 
an early-type population of individually bright galaxies 
with a diminishing faint-end whereas a the late- type pop- 
ulati on is composed of a rising number of faint galaxies 
(e.g.. iFaber et al.l [2007ft. The characteristic luminosity, 
L* therefore depends heavily on the mixture of spectral 
types at any given epoch. Much work has been devoted 
to the K-band LF where the stellar mass-to-light ratio is 
relatively stable and it can thus be used as an indic ator 
for the stellar mass function (e.g., iCole et ail 120011) . It 
is therefore natural to expect to brighten with cos- 
mic time as more mass becomes locked up in low-mass 
stars. In the red/NIR bands, the luminosity evolution 
is typically AM* ~0.5 — 1.0 between redshift 0.1 and 1 
whereas it is much stronger in the UV/blue rest- frames 
indicating higher star formation rates at earlier times. 
Extensive work has been done on the UV LF which is 
largely driven by its importance as a tracer of star for- 
mation rate and assisted by increasing detection rates of 
distant Lyman break galaxies in deep surveys. We see 
My V brighten with increasing redshift and then turning 
over, thus roughly exhibiting the same behavior as the 
derived star formation history (Madau plot). The wide 
redshift range of available UV LF measurements makes 
it the only LF in which a non-monotonic evolution is dis- 
tinctly seen in My V . In all other bands, the evolution of 
the Schecter parameters can be fitted with an analytic 
function to quantify the global evolution, while "wash- 
ing" out outliers in the process. Several authors have 
parameterized the evolution in indi vidual bands (e.g., 
iLin et all 119991 : iCirasuolo et al.ll2010l ). but to our knowl- 
edge, our work is the first multi-wavelength parametric 
study of the evolution of the LF parameters. We find the 
following forms to fit the data well across our wide range 
of wavelengths and redshifts: 

M*(z) = M*(z ) - 2.5 log [(1 + {z- z a )f] (3) 
<t>*(z) = <j>*(z ) exp [-p(z - z )] (4) 

and we assume the following a priori form for the faint- 
end slope 

a(z) = a{z ) (z/z ) r . (5) 

These fits are shown in Figure El For M*(z) and (t>*{z) 
we have taken z =0-8, but zq= 0.01 for a(z). The other 
best-fit parameters are listed in Table [5] Instead of se- 
lecting a preferred LF measurement for a given redshift 
in each band we have chosen to let all measurements con- 
tribute equally to the fitting process regardless of depth, 
area and sample size of the survey. Although there are a 
few notable discrepancies between the data and the fits 
we note that our IR-fluctuation results are unaffected as 
long as the fits remain good in the light shaded areas 
of Figure [3] These regions correspond to the distance 
for which the rest-frame emission is redshifted into the 
observed near-IR wavelengths of interest, defined to en- 
compass the 1.25-4.5/im range. In the following sections 
we will rely on lightcones extrapolated from the highest 
measured redshift, typically z^4, out to z max =7 (see Ta- 
ble [2]). To account for the turnover observed in M*, y , we 



TABLE 2 

Properties of the data shown in Figure \3\ and the best-fit 
evolution parameters of equations ©-([sj 



Band 




N 


%max 




<t>l ,P 


Q0,r 


(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


(7) 


UV 


0.15 


24 


8.0 


-19.62,1.1 


2.43,0.2 


-1.00,0.086 


u 


0.36 


27 


4.5 


-20.20,1.0 


5.46,0.5 


-1.00,0.076 


B 


0.45 


44 


4.5 


-21.35,0.6 


3.41,0.4 


-1.00,0.055 


V 


0.55 


18 


3.6 


-22.13,0.5 


2.42,0.5 


-1.00,0.060 


R 


0.65 


25 


3.0 


-22.40,0.5 


2.25,0.5 


-1.00,0.070 


I 


0.79 


17 


3.0 


-22.80,0.4 


2.05,0.4 


-1.00,0.070 


z 


0.91 


7 


2.9 


-22.86,0.4 


2.55,0.4 


-1.00,0.060 


J 


1.27 


15 


3.2 


-23.04,0.4 


2.21,0.6 


-1.00,0.035 


H 


1.63 


6 


3.2 


-23.41,0.5 


1.91,0.8 


-1.00,0.035 


K 


2.20 


38 


3.8 


-22.97,0.4 


2.74,0.8 


-1.00,0.035 


L 


3.60 


6 


0.7 


-22.40,0.2 


3.29,0.8* 


-1.00,0.035 


M 


4.50 


6 


0.7 


-21.84,0.3 


3.29,0.8* 


-1.00,0.035 



Note. — 1) Fiducial rest-frame band, (2) the effective wave- 
length in microns, (3) number of LFs used, (4) highest redshift 
of LF available in band, (5) Best-fit parameters for M*(z) with 
20=0. 8, (6) Best-fit parameters for <p*(z) with 2o=0.8 in units of 
10 — 3 Mpc — 3 , (7) The parameters for a(z ) ch osen to reflect the 
models (HFE&LFE) presented in Section l3.ll 
'assumed to be the same as in K 



only use our Equation ([3]) out to z^3 w here they inter- 
sect th e high-z fitting formulae given by iBouwens et al.l 
(|2011f ) which we adopt for z>3. 

Evolution is not easily discerned in the faint-end slope, 
a, which by the very nature of surveys is hard to measure 
over large distances. For this reason we explore different 
scenarios for the behavior of a(z) which we explain in 
Section 13.11 In the L and M bands, the redshift range 
covered by the available measurements is so limited that 
we can only fit M*(z) but not the other Schecther pa- 
rameters. Thus, for these two bands we assume <t>*{z) to 
take on the same form as the neighboring A'-band. For- 
tunately, the data available in the LM-bands covers the 
redshift range of interest as is indicated by the shaded 
regions in Figure [3] 

There is significant degeneracy in the Schechter pa- 
rameters derived for a given galaxy population which 
can manifest itself in different values of (a,M* ,<$*) de- 
pending on the LF-estimator used (see Appendix [B]l . 
The overall shape of the LF can appear similar de- 
spite different Schechter parameters typically resulting 
in a comparable value for the luminosity density, C — 
<f)*L*T(a + 2), which we d isplay in the botto m panels 
in Figure IS For exampl e, lllbert et al.l (f2005h (VVDS) 
and lGabasch et all (|2006f ) (FDF) derive comparable LFs 
depite giving very different values for the Schechter pa- 
rameters. The general agreement of the £-data and the 
curves, 4>*(z)L*(z)T(a(z) + 2), indicates that our sepa- 
rate fits do not systematically over- or under-estimate 
the total luminosity density. 

The second step is populating the lightcone seen from 
the standpoint of the observer. Light from distant galax- 
ies appearing in the observed X-hand was emitted at 
wavelength Ax /(l + z) i.e. at all rest- frame wavelengths 
shortwards of Xx throughout the redshift cone. We ex- 
tract the Schechter parameters from our fits in Figure [3] 
at the redshift defined by Zi = Xx/ X Y i — 1 where Y corre- 
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0.1 1 10 0.1 1 10 

z z 

Fig. 4. — Evolution of the Schechtcr parameters and the lumi- 
nosity density seen in the observed at 2.2/im (light squares) and 
3.6/im (dark diamonds). The values are extracted from the fits in 
[3] at the appropriate redshifts. 0* is in units of fO~ 3 Mpc~ 3 and 
vC v in units of 10 40 erg-s — i-Mpc -3 . 



sponds to our fiducial bands (UVUBVRIzJHKLM) and 
Ay < Xx- Our template LFs then become 

*i(M\zi) = 0.4 In (10)0* (zi) (io°- 4 ( M *^)- M )) Q(Z ' )+1 

xexp(-10°- 4 ( M *^)- M )). 

(6) 

The continuous evolution of the LF seen in the X-band 
is then obtained by interpolating the $i's from z = to 
Zmax- It should be noted that because of the a — M* 
degeneracy, our separated (a(z) ,M* (z) ,4>* (z)) fits used 
in Equation ©, cause some amount of deviation from 
the original shape of the LF. This is a small effect in 
comparison with the general disagreement between in- 
dividual authors on the shape of the LF. We refer to 
Appendix [A] where an independent method is used to 
populate the lightcone, in which the original shapes of 
the LFs are kept intact. We show that the two differ- 
ent methods produce the same results, confirming the 
validity of our standard treatment. 

As an example we show in Figure |4] the Schechter pa- 
rameters characterizing the LFs, probing the sky in two 
different observer-frame bands centered at 2.2, 3.6 /xm 
respectively. Although the abundance of galaxies dimin- 
ishes by itself at high-z according to our fits, we impose 
a limit of z max =7 in our modeling, beyond which we 
assume that ordinary galaxy populations were not yet 
established. But due to the steep drop of 4>* at high-z, 
our results are not sensitive to this parameter: in fact, 
using z max — 30, yields results nearly identical to our 
fiducial model. We emphasize that our evolution models 
are empirically supported out to z^4 only, beyond which 
we extrapolate the evolution deduced at lower redshifts. 

In order to deduce the rest-frame LF from survey data, 
absolute magnitudes need to be derived from apparent 
magnitudes. We derive the flux from galaxies in our 
lightcone by backtracking the original procedure i.e. go- 
ing from absolute magnitudes back to apparent magni- 
tudes. This implies undoing any corrections the authors 



have made in this process 

to = M + DM(z) + K{z) + E{z) + A b (l, b), (7) 

where DM(z) is the distance modulus, K(z) is the k- 
correction, E(z) is the evolution correction and At, is the 
correction due to galactic extinction at the Galactic coor- 
dinate (l,b). In LF measurements, authors typically use 
de-reddened magnit udes or correct for ex tinction using 
Galactic dust maps (ISchleeel et al.lli~998l) . This correc- 
tion can be large in the UV/optical but becomes less 
severe towards the infrared where we have Ay /Ak~7-10 
approaching ~15-20 in the IR AC bands (jCardelli et al.l 
1989: (indebetouw et al.l l2005h . We are only concerned 
with emission entering the Milky Way as near-IR where 
the extinction correction is typically well within 0.1 mag 
so we neglect it in Equation (JT]). Correcting for evolu- 
tion is intended to make a sample drawn from a distribu- 
tion of redshifts reflect the true luminosity function at a 
given epoch (usually z me d of the survey/bin) by account- 
ing for ch anges in luminosity a nd number density over 
time (e.g.. lBlanton et al . 2003b). This has been done for 
some local surveys where a considerable spread in the 
rcdshift distribution leaves more cosmic time for evolu- 
tion to take p lace. This typic ally results in corrections 
of ~0.1 mag (jBell et al.l 120031) but since the evolution 
correction simply acts to make the LF more accurate 
at a given redshift we do not need to make any adjust- 
ments. The only magnitude a djustment in Eq uation ([7]) 
of concern is the k-correction (jHogg et al.ll2002l ) which is 
needed to transform to the rest-frame by accounting for 
the redshifted SED of a given source. There are a va- 
riety of methods to deal with this SED dependence and 
we refer to Appendix [B] for a more complete discussion 
of two commonly used treatment in the literature. From 
the k-corrected absolute magnitudes, we simply require 
the spectral independent term to account for the red- 
shift into the observed frame, K(z) = — 2.51og(l + z). 
Equation ([7]) is now reduced to 

m = M + DM(z) -2.5 log(l-M), (8) 

which is the conversion we use. In Section [4] we show 
that we recover the observed number counts to a very 
good accuracy using this methodology 

3.1. The Faint- End LF Regime 

The source subtracted CIB fluctuations are isolated to 
faint sources. By the nature of galaxy surveys, the faint- 
end is generally poorly constrained causing large uncer- 
tainties and scatter in measurements of a, especially at 
high-z. Because of this, many authors prefer to keep a 
fixed in their Schechter fits. Since the data does not show 
robust evolution in a in most bands (unlike M* and </>*) 
we explore variants of the behavior of the faint-end slope 
to get a feel for the sensitivity of CIB fluctuations to the 
abundance of faint galaxies. The substantial scatter in 
measurements of a leaves us some freedom in modify- 
ing the faint-end regime but we find that deep galaxy 
counts impose strict limits on the allowed range of faint- 
end slopes. This is most notable in BVRI, where a steep 
faint-end at z=l-3 leads to an overproduction of the ob- 
served JHK number counts in the faintest magnitude 
bins (see Figure [5]). We therefore consider the range 
of allowed a(z) scenarios that collectively yield galaxy 
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counts consistent with observations across all bands si- 
multaneously. We leave M* and <j>* unchanged when 
varying a despite degeneracies in the parameters (see ap- 
pendix^. We consider two models, high faint-end (HFE) 
and low faint-end (LFE), which, based on the resulting 
galaxy counts, are likely to bracket the true behavior of 
the faint-end of ordinary galaxy populations. These are 
shown in Figure [3] and [5] as the upper and lower bound- 
aries of shaded regions. With the faint-end reasonably 
well constrained at z=0, ranging from -0.8 to -1.2, we fix 
a at these two values for LFE and HFE respectively and 
vary later evolution by changing the slope of the power- 
law in a(z) (called r in Equation ©JE Our HFE model 
is characteristic of str ong steepening such as that found 
by lllbert et all (|2005| ) (VVDS) out to z~l whereas the 
LFE implies a m o re modest e volution, closer to that of 
iMarchesini et all (|2007l l2012j ). Our LFE reflects a lack 
of evolution in the NIR i.e. q~ const., which seems t o 
be favored by some authors (e.g.. iCirasuolo et alj [2007). 
We choose a faint-end cutoff for each template LF at 
L cut = \Q- i L* for LFE and 1CT 8 L* for HFE, thereby ex- 
trapolating the LF to very low luminosities. For both 
scenarios we find 10 -5 L* to be near saturation with 
flux contribution for fainter magnitude bins always being 
<0.02 nWm _2 sr _1 . Our "default" model is the average 
of HFE and LFE with a cutoff at 10 _5 L*. 

We have chosen our LFE/HFE models so that they 
remain consistent with number counts data. The LFs 
dominating the faint counts in Figure [5] are mostly de- 
termined by the faint-end slope, a, at high and inter- 
mediate redshifts and it is important to emphasize that 
more extreme faint-end evolution models generally yield 
number counts that are inconsistent with observations. 
Alternatively, one could in principle imagine an increase 
in the LF in the faintest magnitudes observed deviating 
from a Schechter function. In fact, such an upturn has 
been observed locally, for which a "double" Schechter 
function provides a better overall fit o f the LF (e.g. 
iBlanton et all 120051: lLovedav et all 120121) . Allowing for 
a much steeper slope at z=0 to accommodate this possi- 
bility does not affect the resulting CIB fluctuations be- 
cause the surface density of sources on the sky tends to 
be dominated by populations at larger distances. This 
can be illustrated by examining the underlying LFs of 
the resulting galaxy number counts in Figure [5j where 
the gray lines starting at the bright-end (from left) cor- 
respond to the local contribution (the thick line being 
the most local) moving to high-z LFs to the right. The 
rapid redshift evolution of the cosmic volume element 
prevents a large surface density of low- 2 sources and we 
find the faint counts always being dominated by popula- 
tions at intermediate and high redshifts (z>l). In order 
for low-z sources to have sufficient densities to dominate 
the faint galaxy counts, and thereby also the unresolved 
fluctuations, we would need an extremely steep faint-end 
at z=0, becoming flatter towards increasing redshift i.e. 

igh-z which is the opposite of the observed 
evolution trend. Alternatively, a sudden upturn devi- 
ating from a Schechter form at faint magnitudes would 
need to shoot up by roughly two orders of magnitude. 

1 In the rest-frame UV/optical, where the low-z contribution 
does not matter for the observed NIR, we fix the low-z slope at 
-0.9 and -1.1 for LFE and HFE respectively. 



We therefore consider our HFE scenario to be sufficiently 
extreme at low-z and making it steeper does not have 
an effect on our results. On the other hand, if a sig- 
nificant upturn in the LF exists at z>0.5 (so far unde- 
tected), then this may result in a non- negligible contri- 
bution to the unresolved fluctuations. The large number 
of small halos predicted by the standard ACDM model 
permits such a scenario, especially if the first population 
of dwarfs with normal stellar populations fo r med in halos 
with mass <10 9 M (|Ricotti et al.ll2002allbl . [200cl . For 
instance, if the ultra-faint dwarf galaxies recently discov- 
ered around the Milky Way can be identified as fossils of 
the first galaxies formed before reionization, that would 
imply that we have only discovered a small fraction of 
a widespread population of d warfs which were almost 
certainly brighter in t he pas t ([Ricotti fc Gnedinl 120051 : 
IBovill fc Ricottil 120091 feOllallH ). However, it is unclear 
how to make the flux from this population sufficiently 
large to reproduce the measured fluctuation signal and, 
furthermore satellite dwarfs are efficiently masked along 
with their host galaxy in fluctuation measurements as 
displayed by the masking t ypical ly having angular ra- 
dius of ~15' (lArendt et all i|20Tol ). see also Fig. A-3 in 
IKashlinskv et all (|2012D T In this work we probe whether 
the known galaxy populations, which we extrapolate to 
faint magnitudes in our HFE and LFE limits, can account 
for the observed source-subtracted CIB fluctuations, and 
the question of the nature of the new populations that 
can explain these fluctuations is, while important, out- 
side the scope of the current discussion. 

4. NUMBER COUNTS AND BACKGROUND LIGHT 
FROM LF DATA 

Galaxy number counts have the advantage of being free 
of the uncertainties associated with e.g. k-corrections 
and redshift determination which makes it an important 
test of both cosmology and galaxy evolution models. We 
project our lightcones onto the sky to obtain the galaxy 
number counts in each magnitude bin per unit solid an- 
gle: 

N(m) = J Hm\z)-^dz, ( 9 ) 

where dV/ dzdQ is the comoving volume element per solid 
angle. In Figure [5] we display the number counts from 
Equation (jH]) in the 0.45-4. 5fim range and compare with 
existing data in the literature. The good agreement be- 
tween our modeling and observed counts demonstrates 
the validity of our method. We also display the range 
bracketed by out two limiting models for the faint-end 
slope of the LF, as discussed in Section 13.11 (shaded ar- 
eas) . The gray curves in Figure [5] reflect the underly- 
ing template LFs contributing to the number counts in 
different redshift bins (bright/left to faint/right corre- 
spond roughly to low-z to high-z), elucidating the differ- 
ent populations governing the source surface density on 
the sky. It is reassuring, although not surprising, that we 
recover the shape of the galaxy counts using independent 
observations (the only assumption being the Schechter 
parametrization of the LF). This explicitly confirms that 
our multi- wavelength collection of observed LFs provides 
an accurate description of the photometric properties of 
resolved galaxies on the sky. 
Figure [6] (left) examines how the shape of the number 
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Apparent magnitude (AB) 

Fig. 5. — Galaxy number counts in our default description (solid curve) including the regions of bracketed by our two extreme models, HFE 
and LFE (gray shaded areas). The gray curves show the underlying template LFs in our fiducial bands (Equation which we interpolate 
and integrate to obtain the number counts via Equation The low-z LF dominate the bright counts whereas high- and intermediate 
redshift LFs dominate the faint counts (from left to right). The most local available L F is shown as thick gray curves to demonstrate 
their negligible contributio n to the faint counts. F or 0.45-0.80/xm pa nels the data are from Capak ct al. (2004) (red asterisks) Capa k et al.l 
(2007) (purple diamonds), McCrackcn ct al. (2003) ( green triangles), Yasuda ct al. (2001) (turqoisc diamonds) and|Kashikawa ct al. (2004) 
(blue sq uares). Data in the 1.25-2. 2/rni panels are taken from Vaisancn ct al. (2000) ( green triangles), D ickinson et al. 1 999 (purple 
squares), Maihara ct al. (2001) (green asterisks), Kcenan ct al. (2010b) (blue triangles), Kccnan ct al. (2010a) (blue diamonds) Frith et al. 
(2006) (yellow triangles), Thompson ct al. (2005) (yellow asterisks), Metcalfe et al. (2006) (green diamonds), Quadri ct al. (2007) (turqiose 
triangles), Baker ct al. (2003) (purple squares), Minowa ct al. (2005) (orange squares), Huang ct al. (1997) (red crosses) and the 3.6-4.5/im 
data comes from lFazio et al.l l|2004h (purple symbols). 




Apparent magnitude (AB) Redshift 

Fig. 6. — Left: Our reconstructed number counts in BRIJKHLM compared across the spectrum. The counts have been multiplied by 
a slope of 10 — 4 to bring out features in the shape. In this representation is proportional to the flux contribution from each magnitude 
bin. Right: The accumulation of integrated background light from galaxies over time. The flux builds-up from high-z (right) to low-z (left) 
reaching the present-day observed value listed in Table [3] 
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counts varies across the spectrum (0.45-4.5/xm). Both the 
shape and amplitude of the counts are governed by the 
behavior of the (a(z),M*(z),</>*(,z)) b s -parameters shown 
in Figure |4] and some particular features deserve a few 
remarks. The bright counts all start out with a well 
known (Euclidian) slope of dlog N/dm ~0.6 continuing 
down to to~18-20 where it flattens to ~0.4. To first order 
this "knee" is simply caused by the transition from M*- 
dominated to a-dominated regime. More specifically, a 
dip appears in the BVRIJ number counts at m ~18- 
20 which arises from the lack of very bright galaxies at 
higher redshifts, i.e. M* bs becomes fainter with redshift 
(see Fig. @| . At higher redshifts (and shorter rest- frame 
wavelengths) we see a brightening again which is associ- 
ated with star forming galaxies, bright in UV rest-frames. 
This brightening causes another feature at ~25 mag re- 
vealing a "double-knee" surrounding the dip. This is 
most pronounced in the f?Ui?/-counts but disappears at 
longer observed wavelengths where the UV rest-frame 
becomes too distant. Beyond m~25-26, the counts are 
gradually diminished by the ACDM volume clement. De- 
pending on the exact faint-end model, the logarithmic 
slope in this regime is ^0.2-0.3 in BVRI, decreasing as 
we go to longer wavelengths. 

Another clear feature of the number counts seen in 
Figure H] (left), is the overall increase per magnitude bin 
as we go to longer observed wavelengths. We find the 
reasons for this to be twofold. First, the bright end is 
typically dominated by galaxies which are more lumi- 
nous in the red bands such as the case of giant ellipti- 
cals. Therefore we see a larger number of them out to 
greater distances (in Fig. [3] we clearly see M* becoming 
overall brighter from blue to red). Second, when we look 
at the Universe through redder bands, we observe the 
redshifted light from bluer rest-frames emitted in epochs 
when the star formation activity was greater and conse- 
quently M* was brighter. We further point out that our 
reconstructed counts are immune to confusion and agree 
w ell with the confu sion corrected Spitzer/IRAC counts 
of lFazio et al.l (|2004[) (confusion enters around toab~20- 
22). 

We infer the amount of background light from galaxies 
from our reconstructed counts: 



vL, 



f(m)^-dm, 
am 



(10) 



where /(to) = vf v of Equation @ and T is the in- 
tegrated flux in units of nWm _2 sr _1 . Figure [6] (right) 
shows how extragalactic background light builds up with 
cosmic time observed through BVRIJKLM. This re- 
sults in present day values of the integrated background 
light of 9.6, 9.3, 8.1, 4.9 and 3.3 nWrn^^r" 1 at 1.25, 
1.63, 2.2, 3.6 and 4.5^m respe ctively (see TableED ), which 
agree very well with Table 5 of Kashlinskv (2005) and are 
also in general agreement with lMadau fc Pozzettil (12 000). 
but slig htly lower than the values found by lKeenan et all 
(|2010aT ). A subtle underestimation could be due to the 
smooth fitting of the LF evolution which smears out 
any abrupt variation of the Schechter parameters which 
could either be physical. The small d eficit with respect 
to the EBL of IKeenan et al.l (|2010aj) arises in the 21- 
23 mag range wh e re we see a better agreement with 
iMadau fc Pozzettil (|2000f ) and lMaihara et all (|2001h . 



5. NEAR-IR FLUCTUATIONS FROM 
UNRESOLVED GALAXIES 

We now turn to evaluating the source-subtracted CIB 
fluctuations keeping in mind the procedure leading to 
their detection from raw images. If enough pixels re- 
main in the maps after the masking of resolved sources, 
the fluctuations can be characterized via their angular 
power spectrum, which can then be computed more effi- 
ciently by using FFTs than the 2-point correlation func- 
tion. For a detailed description of the process of reducing 
CIB fl uctuation dat a in th e Spitzer/IRAC analysis we re- 
fer to lArendFeFal] (11013) . 

The measured two-dimensional power spectrum from 
extragalactic sources consists of two components: i) 
the shot noise from the fluctuation in the number of 
unresolved sources entering the instrument beam, and 
ii) the clustering component arising from the correla- 
tion of galaxies on all scales. Additional power arising 
from local components such as Galactic cirrus and Zo- 
diacal Light has been shown t o be comfortably below 
the measured signal a t 1-5/xm (Kashlin skv et al.l 120051 : 
iMatsumoto et al.l[20TTI : iKashlinskv et al.ll2012h . In com- 
paring with observational data we adopt the conven- 
tion for the power spectrum to approximate the root- 
mean-square fluctu ations as (q 2 P 2 (q)/2Tr) 1 / 2 ~ (5F 2 ) 1 / 2 
(|Kashlinskvl I2005D . The angular power spectrum of 
galaxies projected onto the sky can be related to their 
evolving 3D power spectrum, P^(k), by the Limber ap- 
proxim ation (for 9 <Cl radian) which we adopt as modi- 
fied bv lFernandez et al.l (|2010l) . 

1 f \dF~\ 2 Ps(qd^ L 1 ;z) dz 



dz 



MM 



1 + z' 



(11) 



where dA is the comoving angular diameter distance. 
The quantity in the square brackets is the flux production 
rate which is empirically determined by our populated 
lightcones: 



dT 

dz 



dmf(m) 



dN{m\z) 
dz 



(12) 



It is important to note that the process developed in 
IKashlinskv eTaTl (120051 [2007b) removes sources down to 
a fixed level of the shot- noise power (see Table HJ. This 
is equivalent to removing galaxies down to a limiting 
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All quantities are in nWm~ 2 sr~ 1 . 
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Fig. 7. — Flux production rate (times z) as a function of red- 
shift in the unresolved regime shown for limiting magnitudes of 22, 
24, 26, 28 (solid, dashed, dot-dashed, dotted curves respectively). 
The total unresolved flux under each curve listed in Table E] The 
figure illustrates how removal of ever fainter sources isolates the 
unresolved component to higher redshifts. 

magnitude, mu m , so that the remaining unresolved back- 
ground is given by Equation (|12[) integrated from mi im to 
oo. In Figure [7j we show the unresolved background from 
our modeling as a function of redshift, which illustrates 
the process of galaxy removal down to fainter magni- 
tudes isolating the background to progressively higher 
redshifts. Note, that there is very little contribution 
(<0.1 nWm _2 sr _1 ) from galaxies at z<l after removing 
galaxies down to 26 AB mag. We find that for a limiting 
magnitude brighter than ~24 mag, the unresolved flux 
is mostly dominated by M* galaxies at intermediate red- 
shifts whereas galaxies at the faint-end takes over once 
mii m >24:. In Table [3] we list the total integrated back- 
ground in the 0.45-4.5^m range including the unresolved 
background for different limiting magnitudes correspond- 
ing to the curves in Figure [71 

5.1. Shot Noise 

The shot-noise level seen in fluctuation measurements 
is critically important in o rder to identify the natur e of 
the unresolved populations Kashlins kv et al.l (|2007d) . It 
can be described as statistical counting noise in the num- 
ber of unresolved sources within the instrument beam 
and its power is, 



SN 



dz 



dm f 2 (m) 



dN{m\z) 
dz 



(13) 



Shot noise is a directly measurable quantity and is not 
affected by confusion which may be present. This allows 
us to evaluate the effective limiting magnitude, mi im , for 
a given shot noise level using our models which are also 
immune to confusion. We calculate the shot noise associ- 
ated with galaxies in our lightcones and display it in Fig- 
ure [5] as a function of limiting magnitude at the relevant 
bands. As fainter galaxies are removed the shot noise 
drops steadily in the same manner as seen in measure- 
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Note. — The upper and lower values are not error but corre- 
spond to the HFE/LFE evolution scenar ios of the faint-end slope . 
*The values are inferred from Figure 3 of [Matsumoto ct al. 120111 ). 

ments. At ^22 mag we have already removed most M* 
galaxies at z<l beyond which the shot noise is mostly 
determined by the faint-end of the LF. The horizontal 
lines in Figure [5] show the levels reached by the stud- 
ies listed in Table |H The intersection with our mod- 



els agrees well with [Kashlinskv ct al. (2007b) claiming 
to have removed galaxies down to m ^25-26 AB mag 
bu t is slightly brigh t er (m ~24) for the levels reached 
by iKashlinskv et al] (|2005l ) who claimed to reach ^25 
mag. Si milarly, our shot noise leve ls agree well with those 
found by iMatsumoto et all ([20 111) after removing galax- 
ies down to AB magnitudes 22.9, 23.2 and 23.8 in the 
AKARI/IRC bands at 2.4, 3.2 and 4.1^m respectively. 
Table [4] lists the limiting magnitude predicted for the 
shot noise levels reached in several studies. 

We have defined mu m to separate resolved/removed 
galaxies from unresolved remaining sources. In practice, 
however, the accurate value of mu m reached depends on 
the source detection algorithm and the photometric aper- 
ture used to derive magnitudes. Furthermore, source ex- 
traction can become limited by confusion, depending on 
exposure and instrument beam. Since our underlying 
reconstruction of galaxy counts from LFs is immune to 
confusion, we assume that the measured shot noise lev- 
els serve as a reliable indicator for the faintest sources 
removed, mu m . This obviously assumes that the source 
removal is done properly and does not introduce spuri- 
ous signa ls in the backg r ound fluctuations as discussed at 
length in lArendt et al.l ([2010D . It also assumes that the 
(quasi-)flat power seen on small scales is entirely due to 
shot noise dominating the contribution from non-linear 
clustering of galaxies which we discuss in the following 
subsection. 

5.2. Galaxy Clustering 

The shape and amplitude of the fluctuations produced 
in each redshift slice is dictated by the two evolving 
quantities in the Limber equation (eqn. 1111) . i) the 
amount of light production given by our reconstructed 
dJ- / dz in a given band, and ii) the clustering pattern 
of the sources in this epoch, described by their three- 
dimensional power spectrum, P^{k,z). For the latter 
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Fig. 8. — Shot-noise power amplitude after integrating the counts as a function of limiting magnitude (connected squares). The gray 
shaded area corre s ponds to the all owed range of the faint-e nd evolutio n of the LF. The thick gr ay lines show the levels of P$n reached by 
IMatsumoto et al.l fed 111) at 2.4um. lKashlinskv et al.l (120051) (dark) and lKashlinskv et all {2007b) (light) at 3.6 and 4.5/im. The intersection 
corresponds to the limiting magnitude reached in these studies. We tabulate these values in Table [3] We point out that our model counts 
are immune to the effects of confusion. 



quantity we assume that on large scales sources clus- 
ter according to the observationally established con- 
cordance ACDM power spectrum. Prescriptions exist 
for non-linear evolution that modify the linear power 
spectrum in the regime where structures have collapsed 
out of the density field and linear theory b reaks down 
(peacock fc Doddsl Il996l ; iSmith et al.l l200l . However, 
luminous sources are known to be biased tracers of the 
dark matter distribution particularly in the non-linear 
regime where the correlations of sources depends on the 
Halo Occupation Distribution (HOD) of galaxies. We 
therefore consider a halo model description of the power 
spectrum which decomposes it into two terms, a two- 
halo term (P 2h ) on large scales arising from the cor- 
relations of isolated halos, and a one- halo (P lh ) from 
correlations of particl e s with in the same halo on small 
sc ales (e.g. iMa k Frvl ti2000D '). We follow the treatment 
of iCoorav & Shethl (|2002t ) and write, 



We take the mass dependence of ou r HOD mode l to fo l- 
low the four parameter description of Zheng et al. ( 2005): 



(N cen ) 



(Nsat) = J 



1 + erf 



1 + erf 



logAf -logAf mj 

ClogM 

log M - log 2M„ 

ClogM 



(18) 



M 



(19) 



P 9al (k) = P Ul (k)+P 2h (k), 



(14) 



where, 



P lh (k) 



dM 



dn 2(N sat )(N cen )u(k\M) + (N sat ) 2 u 2 (k\M) 



dM 



gal 



P Zh (k) = P tm (k) 



dM 



dn (Ng a i) 



dM 



ngal 



b{M)u{k\M) 



(15) 



(16) 



and dn/dM is the halo mass function (jSheth et al.ll200il 
from), fig a i is the average number density of galaxies, 
P Un (k) is the linear ACDM po wer spectrum (c o mput ed 
using the transfer function of Bardccn et al.l ([1986)), 
u(k\M) is the normalized F ourier transform of the halo 
profile dNavarro et al.|[l9 96V and b(M) is the halo bias 
(from ISheth et al.l 1200 lf T The occupation number has 
been separated into central galaxies, (N cen ), and satel- 
lite galaxies, (N sat ), such that 



(Ngal) = (N cen ) + (N sat ). 



(17) 



where (N cen ) is characterized by M m i n , the minimum 
halo mass that can host a central galaxy and o~\ os m, 
which controls the width of the transition of the step 
from zero to one central galaxy. The satellite term has 
a cut-off mass which is twice as large as the one for cen- 
tral galaxies and grows as a power-law with a slope of 
a s , normalized by M sat . This form has been explored 
both numerically and observationally. Since the mea- 
surements of HOD-parameters are obtained from sam- 
ples of resolved galaxies at low-z, their validity may not 
extend into the unresolved regime or, in particular, to 
higher redshifts. Since we are concerned with the unre- 
solved regime it is important to note that the measured 
cut-off mass of central galaxies, M m i n , is typically set 
by the lowest luminosity probed by the survey so ha- 
los may continue to host central galaxies to lower masses 
but are excluded due to selection criteria. In Section|3]we 
showed how the unresolved light is typically dominated 
by the faint-end of the LF for m>25 with most bright 
central galaxies removed out to z^3 in measurements of 
CIB fluctuations. One would also expect the masking 
to eliminate most of the surrounding satellite galaxies. 
We have adopted the following parameters of the HOD- 
model motivated by SDSS measurements of lZehavi et al.l 
<203): o- losM = 0.2, M min = lO 9 M , M sat = 5 • 1O 1O M , 
and a s = 1 where we have deliberately chosen a lower 
cut-off reflecting low mass halos hosting galaxies well 
into the unresolved regime, and a lower M sat allowing 
for large amounts of unresolved satellite galaxies, while 
keeping a s — 1 . It should be noted, that in the absence 
of any HOD-assumptions, a simple linear ACDM cluster- 
ing with typical bias, b 2 P hn (k), produces nearly identi- 
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cal fluctuations on large scales. The one-halo term has 
white-noise power spectrum (P =const) with its ampli- 
tude limited from above by the measurements at small 
scales and so its modeling is irrelevant to interpreting the 
clustering signal at scales > 1'. 

We assume that unresolved sources in our lightcones 
are uniformly mapped onto the halo distribution i.e. the 
clustering is independent of galaxy luminosity. In prac- 
tice however, we expect the most luminous galaxies to 
be removed in the masking process along with most of 
the accompanying satellites. This could motivate one to 
introduce an upper mass limit in the integrals in Equa- 
tions (fT5|) , M max (z). However, this would require an 
additional mass-to-light ratio assumption and since it 
would always result in a decrease of the clustering ampli- 
tude, we do not apply M max (z) and consider the result 
to be an upper limit for the resultant power spectrum. 
This includes the mass-dependent bias which is simi- 
larly integrated over the entire range of occupied halos 
(^, 10 9 )- The large sca le (linear regime) galaxy bias seen 
bv lZehavi et all ([2011D in the local SDSS sample is b «1 
when all galaxies are includ ed. At somewhat higher red- 
shifts, [GranetFeFal] HoU) find b = 1.38 ±0.05 averaged 
over 0.5<z<1.2. Further increase of the linear bias with 
rcdshift is expected on theoretical grounds as collapsing 
density peaks were increasingly rare in the past. The 
bias pr escription used he re shows the same general be- 
havior (jSheth et al.ll200lD . Several CIB studies at far-IR 
wavelengths clai m a linear bias as high as &=2-3 for far- 
IR sources (e.g.. ILagache et all [2007E IViero et alll2009f) 
but at these redshifts, the samples are already biased to- 
wards the most luminous objects due to selection effects. 
If anything, we expect the bias to be lower in the faint 
and unresolved regime after the more strongly biased lu- 
minous galaxies are masked and removed. 

The large scale fluctuations are always dominated by 
clustering in the linear regime (two-halo term). On the 
other hand, the non-linear clustering described by the 
one-halo term in Equation (ITBl exhibits a P(fc)=const 
behavior (SFocq) making it indistinguishable from shot 
noise in measurements. Given that we found excellent 
agreement between the shot noise in our models and 
the measurements at the same magnitude levels, there 
does not seem to be any need to invoke non-linear clus- 
tering to explain fluctuations on small-scales (unless the 
data points deviate from a simple white noise spectrum 
SFoaq). In addition, we explored the pure dark-matter 
treatm ent of the non-linear clustering of iSmith et al.l 
(2003) but find it to be inconsequential in compari- 
son with the shot noise dominated fluctuations on small 
scales. Although we see the one-halo term contributing 
somewhat to the HFE fluctuations in FigureEE it becomes 
less relevant if one accounts for the more massive halos 
being masked/removed. In fact, we will see in Section[5]4] 
that current fluctuation measurements place a limit on 
the amount of non-linear power in the unresolved regime. 

5.3. Comparison with fluctuations from the 
Millennium Simulation and Semi-Analytic Models 

To compare our results with the clustering of ha- 
los seen in large scale N-body simulations, we have 
ma de use of the theore tical lightcones constructed 
by iHenriques et al.l (|2012l ). These mock catalogs are 
based on semi-analytical models for galaxy evolution 



(|Guo et al.l 120111) which are implemented on two very 
large dark matter simulat ions, the Millennium Simula- 
tion dSpringel et al.l I2005T) and the Millennium-II Sim- 
ulation (jBovlan -Kolchi n et al.l [2009). The simulations 
provide a description of the evolving spatial distribu- 
tion of dark matter halos and subhalos whereas the na- 
ture of the baryonic content is described by the latest 
versio n of the semi-analytical Munich model (|Guo et al.l 
120111 ). The Millennium Simulation follows structure for- 
mation in a box of side 500ft.~ 1 Mpc comoving with a res- 
olution limit of ^lO lo /i _1 M whereas the Millennium- 
II Simulation focuses on a region of 100/i _1 Mpc but 
with complete merge r trees down to ^lO 8 /i _1 M G 0. 
IHenriques et al . (2012) use the Millennium Simulation 
only in their study, limiting the faint-end of the LF 
to halos >lO lo /i _1 M0. Even so, the predicted faint 
near-IR counts are higher than observations suggest due 
to an unusually high abundan ce of relati v ely lo w mass 
galaxies (~1O 1O M ) at z>l (|Guo et al.l (|2011D tuned 
their model to match the local populations). A com- 
parison of the predicted correlation function of these 
models with local SDSS data shows decent agreement 
for massive galaxies whereas correlations of low mass 
syste ms are overpredic t ed, pa rticularly at small separa- 
tions. IHenriques et al.l (|2012D also neglect the effects of 
dust and PAH emission and consider only starlight us- 
ing s tellar synt h esis m odels of IBruzual fc CharlotJ (|2003D 
and iMarastonl (|2005[ ). For a detailed de s cription of 
these mode ls we refer tolSpringel et all (|2005l) , IGuo et ah! 
(|20T1 and IHenriques et all (|201^T 

Despite the limitations mentioned above, we find that 
this study provides a useful comparison to our fluctuation 
analysis. After constr ucting images using the publicly 
available mock data of Hcuriqu es et al.l (|2012l ). we cal- 
culate the projected angular power spectrum, convolved 
with the instrument beam. We analyze two independent 
regions observed in H, K, IRAC1 and IRAC2 each cover- 
ing 1.4x1.4 degrees on the sky. We extract all galaxies in 
the magnitude range mu m <m<30 to produce the unre- 
solved fluctuations which we display alongside our results 
in Figure[Hl Because of the overabundance of faint galax- 
ies at 3.6 and 4.5 jLtm in the semi-analytical description of 
IGuo et alJ (f20"TTh . we need to remove galaxies down to 0.2 
mag deeper than the mu m listed in the panels in order to 
normalize to a common shot noise level. This NIR over- 
abundance (despite the resolution limit of ^10 10 h~ 1 M & ) 
results in the Millennium fluctuations (dark-gray shades 
in Figure [9]) being in closer agreement with our HFE sce- 
nario at 3.6 and 4.5/zm but can otherwise be considered 
to be consistent with our main results. 

5.4. Results 

With the emission history reconstructed from LFs and 
the sou rces distributed according to the halo model in 
Section I5~2l we projected onto the sky the clustering pat- 
tern in our NIR lightcones using Equation (fTTj) and dis- 
play the results in Figure [5] The limiting magnitudes 
have been chosen such as to normalize the shot noise 

2 The Millennium Sim ulation and the resulting lightcones 
of IHenriqu es et al .l H2012T ) assume a WM API-based cosmol- 
ogy j f^pergeT^^^L r i2003l ) with parameters h= 0.73, f2 m =0.25, 
Qa=0.75, n=l and <rg=0.9 which are slightly different that our 
adopted parameters of h= 0.7, f! m =0.3, Qj^=0.7 but this is of no 
appreciable consequence for the results in Figure [9] 
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Fig. 9. — Models of the unresolved near-IR fluctuations compared to measurements from authors listed in the panels. We have chosen 
the limiting magnitude such that the models are normalized to the shot noise levels reached in these studies (including a contribution 
from a one-halo term) . The solid curves show the total contribution from clustering and shot noise whereas the light shaded areas indicate 
the region bracketed by our HFE and LFE models. These are all suppressed by the instrument beam on small scales. The dotted lines 
indicate the separate one-halo and two-halo terms of the power spectrum. Shown in each panel is the total unresolved flux associated with 
the default model (J 7 ), the values of Psn (in units o f nW 2 m sr —1 ) and th e associated m; im . The dark shaded regions correspond to 
fluctuations arising from galaxies in the lightcones of Hcnriqucs et al. (2012) derived from the Millennium Simulation in the magnitude 
range mn m <m<30. Because of their overabundance of faint galaxies at 3.6 and A.5fim we have increased the m; im of t he Millennium 
fluctu ations by 0.2 mag to normalize to the correct shot noise levels. In the 3.6/^m panel we also show the default model from lSullivan et al.l 
(2007) (dashed line). In the 1.6/Ltrn panel the notation follows Fig. 2 of Thompson et al. (2007a): asterisks correspond to fluctuations with 
all sources removed whereas the triangles indicate their estimate of the instrumental Gaussian noise. 
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Fig. 10. — The lines show the ratio of the m easured source-subtract ed power spectrum from the AKARI 2.4/xm data and the latest 
Spitzer-based measurements at 3.6 and 4.5/im (Kashlinskv ct al. 20l2) to the HFE and LFE expectations (red/upper and blue/lower 
respectively). The results show that the measured CIB fluctuations continue to diverge from our models as we go to larger scales and are 
thus unlikely to result from extra-biasing of these faint populations: in order to explain the measured signal the biasing would have to be 1) 
scale-dependent, i.e. non-linear, 2) biasing amplification would have to be more non-linear on scales where the amplitude of the underlying 
correlation function is weaker (larger scales), and 3) the biasing would have to be different at 3.6 and 4.5/tm. 



(dot-dashed lines) to the measurements shown in each 
band. The shot-noise is seen to dominate the fluctua- 
tions on small scales whereas the clustering component 
becomes significant at arcminute scales. In the display we 
have chosen to focus on 1.6, 2.4, 3.6 and 4.5/xm where we 
can compare with measurements from Hubble /NICMOS, 
AKARI/IRC, and Spitzer/IRAC. Our models have been 
convolved with the beam profile (or PSF) of these instru- 
ments. It is immediately clear from Figure 03 that the 
contribution from known galaxy populations falls short 
of the measured clustering signal in every band shown. 
We briefly discuss each com parison: 

IKashlinskv etaLl (I2007bfl find excess fluctuations of 
5F~0. 05-0.1 nWnr 2 sr _1 at arcminute scales in the 
Spitzer/IRAC channels after removing sources down to 
~25 mag or shot-noise levels Psn ^ 3 x 10~ n nW 2 /m 4 /sr. 
It can be seen from Figure [5] that the known sources 
remaining at the measured shot-noise levels cannot ac- 
count for the observed fluctuations for any faint-end 
modeling of the LF. W e ha ve displayed the data o f 
IKashlinskv eTall (f2005f ) and IKashlinskv eTall (|2007b[ ) 
in panels side-by-side illustrating that the discrepancy 
gets larger as galaxies are removed to deeper levels. 
The unresolved flux associated with our default model 
is 0.18 nWm~ 2 sr~ 1 in the deepest 3.6/im maps of 
IKashlinskv et all (|2007b| ). so in order to explain the ob- 
served level of the excess fluctuations the relative levels 
of the source-subtracted CIB fluctuations would have to 
be close to non-linear, 6F/F~1, all the way to ~ 10'. The 
spatial spectra of the CIB fluctuations from the known 
galaxy populations is such that the gap increases toward 
large scales if this behavior of the source-subtracted CIB 
fluctuations continues as observed (say, ~ 1°), so these 
fluctuations would have to be in the same (quasi)non- 
linear regime at much larger scales making it more diffi- 
cult to explain them with the known galaxies. The addi- 
tional linear biasing to amplify the arcminute scale signal 
to the observed levels but this would require b ~ 6 — 20 
which is highly unlikely for small systems in the 1 < z < 3 
range where most of the flux is produced. This, how- 
ever, can be shown not to be viable in light of the 



latest Spizter-based res ults submitted after this paper. 
IKashlinskv et al.l (|2012t ) measure the source-subtracted 
CIB fluctuations on sub-degree angular scales confirming 
the earlier results and identifying, for the first time, the 
fluctuations spectrum to ~ 1° where the discrepancy con- 
tinues to grow. Figure [TU] shows the ratio of the measured 
power s pectrum from the new l arge scale Spitzer/RiAC 
data of IKashlinskv et al.l (|2012t ) to the power spectra of 
our HFE and LFE (red and blue), illustrating that the 
data keeps diverging from our models out to ~ 0.5°. 
This shows that if one were to model the measured CIB 
fluctuations with extra biasing of the known galaxy pop- 
ulations, the biasing would have to be 1) highly scale- 
dependent, i.e. more prominent on larger scales, where 
density pattern is in linear, 2) the resultant biasing fac- 
tors would have to be huge reaching amplifications of 
over two order of magnitude at the largest scales, and 3) 
the biasing would have to be wavelength dependent at- 
testing to the different discrepancy ratios in each band. 
This argues further against the detected CIB fluctuations 
arising from the a faint-end extension of the known pop- 
ulations. 

Following the HOD- model in Section I5T21 we find that 
small scale power from the one- halo term (P lh ) is lower 
than the shot noise (Psn)- For the HFE model however, 
the two are comparable at 3.6 and 4.5/im which leads to a 
slight increase in the upper values of mu m in Table [4] (by 
0.3 mag) as normalized by the measured shot noise levels. 
The fact that our shot noise, calculated from number 
counts, i s in good agreem e nt wi t h the small scale data 
points of Kashlinskv eTall (p005H2007bh . argues against 
a significant clustering on small scales. It must therefore 
remain at or below the measured level of Ps n > otherwise 
sources would have to be removed to mu m >26, which is 
n ot possible in the de epest Spitzer/IRAC maps. 

iCharv et al.l ()2008l ) stack deep Spitzer exposures to 
detect faint ACS galaxies beyond the dete c tion th resh- 
old of the frames used in IKashlinskv ct al. (2007b]) and 
explore the sensitivity of the IR-fluctuations to these 
ACS sources. Their stacked source detections down to 
26.0-26.2 mag imply a net flux of 0.12-0.35 nWm~ 2 sr _1 . 
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Fig. 11. — The contribution of different redshift bins to the unresolved IR-fiuctuations shown in Figurc[9]for the mu m indicated in the 
panels. The 3.6 and 4.5/xm panels correspond to the models at the shot noise levels of Kashlinskv et al. (2007b). The different set of lines 
correspond to the redshift bins indicated in the legend. This illustrates that depending on the observed band and the depth of source 
removal, the unresolved fluctuations from known galaxies are dominated by populations at different epochs. The amplitude and shape is 
governed by 1) the flux production history (see Fig. [3, and 2) the evolving power spectrum, Ps(k,z). The non-linear clustering component 
is imp orta nt at low- z but moves towards small scales for higher z. The dependence on the comoving angular diameter distance, dj\(z) (see 
Eqn. 111! ) is easily seen as the peak of the ACDM power spectrum shifts towards smaller scales with increasing redshift. 



For comparison, the flux associated with our lightcones 
in the 25-26.2 mag range is 0.04 and 0.2 nWm _2 sr _1 
at 3.6/zm for LFE and HFE respectively with 0.04-0.35 
nWm~ 2 sr ~ 1 from still faint e r galax ies, >26.2 mag. We 
note that iKashlins kv et all (|2007af) demonstrate obser- 
vationally the negligible correlations on arcminute scales 
between the source-subtracted CIB maps, as cons t ructed 
by their self-calibration procedure lArendt et al.l (|2010f ) 
a nd ACS source m aps. 

iThompson etail (|2007aD measure fluctuations at 
1.6/xm on sc ales out to 80" u sing a ST/NICMOS (and 
at 1.1 /jm in IThompson et ahl (|2007bl )) and ascribe the 
signal to faint galaxies emitting at redshifts z~0.5- 
1.5. Their fluctuations at 80" have amplitudes of ^0.4 
nWm _2 sr~ 1 , which is a factor of 2-7 times higher than 
the total unresolved component, 0.06-0.20 nWm _2 sr _1 , 
for sources fainter than >28 mag, indicating that the 
clustering of the underlying galaxies must be highly non- 
linear. For their CIB fluctuation levels to be reconciled 
with our empirical estimates, the one-halo term would 
have to be significantly higher, but then its amplitude 
would overshoot the data at all the other NIR wave- 
lengths. If we take the upper limit on the shot noise 
at these wavelength s to be at the leve l s of the estimated 
instrument noise of IThompson et aTl (|2007a| ). then our 
shot noise already matches at AB magnitude of ~27 (see 
triangles in Figure [S]). But even at that level we cannot 
reproduce the fluctuations (asterisks) with the clustering 
of known galaxy populations out to 1'. We point out in 
this context the clearly visible outer halos of the sources 
removed by Thompson et al (2007a, see their Fig. 4) 
whose contribution to their CIB fluctuations shown may 
be significant and should be estimated for more quanti- 
tative conclusions at 1.6 um. 

The Matsum oto et all (j 2 1 lh measured fluctuations at 

3 We note that in the context of [Thompson et al. (2007aj), our 
theoretical magnitude limit, mn m at 1.6/um should no longer be 
taken as a definitive boundary between resolved and unresolved 
sources because ACS images at shorter wavelengths were used to 
remove sources which can translate to a wider spread in magni- 
tudes at 1.6^tm (due to different exposures and different SEDs of 
individual sources). 



2.4, 3.2, 4.1/im using data AKARI satellite and conclude 
that they are consistent with stars from early epochs con- 
firmin g the identification proposed in [Kashlinsk v et all 
(2005). The left panel in Fig. [TU] confirms that the 
AKARI signal at 2.4/im cannot be explained by the re- 
maining known galaxy populations. 

In Figure [9] we als o display the default model from 
iSullivan et al.1 (j2007P ) (dashed lines) who combined a 
halo model and conditional luminosity functions to cal- 
culate IR-fluctuations at 3.6/im. Our models have a 
somewhat lower amplitude considering the fact that we 
use mu. m =24A as o pposed to the 2 5.3 mag used by 
ISullivan et al.l (|2007l ) (and quoted in IKashlinsk v et al.l 
(200|)) but the two are in rough agreement. For 
mnrn = 25.3 our unresolved flux is 0.1 nWm _2 sr _1 (LFE) 
which is roughly consi s tent w ith the 0.08 nWm _2 sr _1 
found bv lSullivan et al.l (120071). However, they claim t hat 
the fluctuations measured by IKashlinskv et al.l (|2005l ) at 
3.6/im can be explained by galaxies in the magnitude 
range 25.3 to 28.8 (AB) at z^l-3. This is a somewhat 
puzzling conclusion when comparing their model with 
the data in Figure [9] as it clearly fails to account for the 
clustering excesfl 

In Figure [TT] we show the contribution of different red- 
shift bins to the unresolved IR-fluctuations for the mu m 
indicated in the panels of Figure GO This illusrates the 
different epochs in which unresolved galaxy populations 
contribute to the fluctuations in different observed NIR 
bands. The redshift dependence is governed by 1) the 
flux production history (see Fig. [7]), and 2) the evolving 
power spectrum, P^(k,z). The Figure also reflects the 
dependence on the comoving angular diameter distance, 
(Ia{z) (see Eqn. (fTTjl ) with the overall clustering pattern 
shifting towards smaller scales with increasing redshift. 

6. SUMMARY AND DISCUSSION 

We have reconstructed the emission histories seen in 
the near-IR of present-day observers to model the unre- 
solved CIB fluctuations and compared with current mea- 

4 The data-points from Kashlinskv et al. (2005) appear only in 
the electronic version of [Sullivan et al. ( 2 0071) . 
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surements. Our compilation of 233 luminosity functions 
used to populate lightcones at z<7 reproduces the ob- 
served number counts remarkably well and accounts for 
the features shaping them. We assume the Schechter- 
type LF and model the evolution of its parameters from 
the available datasets. We then considered high and low 
faint-end LF limits within the constraints permitted by 
deep galaxy counts data. Extending these to faint mag- 
nitudes and to high-z we calculated the range of unre- 
solved background flux in deep images and derived CIB- 
fluctuations from these galaxy populations predicted by 
the standard ACDM clustering power spectrum. We find 
good agreement between the predictions of our analysis 
and semi-analytical galaxy evolution models combined 
with the large scale Millennium N-body simulation. 

By varying the limiting magnitude of source subtrac- 
tion we normalize our models to the observed shot noise 
levels, finding good agreement with the depths reached 
in current fluctuation measurements. We show that the 
known galaxy populations fail to account for the observed 
source subtracted CIB clustering signal in either LFE or 
HFE limits. Although, in principle, by varying mii m one 
can find a population of brighter galaxies that matches 
the measured clustering amplitude at some fiducial an- 
gular scale, the associated shot noise levels always im- 
ply that all such populations have been removed in the 
source subtraction thereby not contributing to the un- 
resolved fluctuations. Thus it means that the emitters 
producing the source-subtracted CIB fluctuations on ar- 
cminute scales are below the detection limits of current 
surveys and furthermore, cannot be a part of the known 
evolving galaxy populations. In other words, the only 



way to reproduce the clustering excess with extragalac- 
tic sources is by introducing a new population of sources 
that are significantly fainter than the detection threshold 
of current instruments i.e., a highly clustered population 
with low shot noise. 

The high isotropy of the CIB fluctuation s i gnal m ea- 
sured in Spitzer IRAC data Kashlin skv et al.l (|2012T ) ar- 
gues strongly against the signal originating in Galactic or 
Solar system foreground emissions as well as very local 
extragalactic sources. Since the observed galaxy popula- 
tions (extrapolated to very faint limits) cannot explain 
the measurements, the CIB fluctuations must originate 
in new populatio n s so fa r unobserved in galaxy surveys. 
iKashlinskv et all ()2007af) show that there are no corre- 
lations between the ACS maps with sources down to AB 
mag of c^28, and the sou rce-subtracted CIB maps from 
IKashlinskv et al.1 (|2007bf ). This implies that either the 
CIB fluctuations originate in a large unknown popula- 
tion of very small systems at low/intcrmcdiatc rcdshifts, 
or they are produced by high redshift, z>7, populations 
whose Lyman break (at rest 0.12/im) is shifted passed 
the longest ACS channel (at 0.9/xm). 
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APPENDIX 

A. LF BINNING AND INTERPOLATIONS 

Because of degeneracy in (a, M*, 4>*), different sets of Schechter parameters can represent LFs of very similar shapes. 
The method used in Section [3] disentangles the Schechter parameters to fit their evolution individually. In addition to 
this, we used an alternative approach in which the shape of each measured LF is kept intact. We took each LF in its 
rest-frame and redshift the associated emission to the observed wavelength, A obs = A res *(l + z). We examine the all 
LFs that meet the criterion Ao— AA < \ obs < Ao+AA where Ao is the center of the NIR band and AA is roughly the 
FWHM of the filter. The inserts in Figure 10 show the redshift distribution of available LFs which can be observed 
through JHKL. In a given band, we place each LF in redshift bins and take the functional average of <!>(M) in common 
bins so that we have a unique LF at each redshift. We thus have template LFs, $i(M\zi), in each of the observed 
NIR bands and the rest of the analysis is identical to that in Section [3] following from Equation (J6j> (we interpolate the 
evolution and project the populations onto the sky). The major shortcoming of this method is the redshift information. 
Averaging over several LF in a common redshift bins is immune to the effects of Schechter parametrization but comes 
at the cost of crude evolution i.e. the sampling of z is determined by the number of z-bins. As seen in Figure [Al there 
is no guarantee that there exists a LF measurement falling into Ao— AA < A obs < Ao+AA in each redshift bin. In this 
case we borrow LFs from neighboring wavelengths scaling them according to synthetic spectra. Figure 10 shows that 
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Fig. 12. — Comparison between our default method (dashed) and the alternative method presented here (solid). The two curves agree 
to within 20% in the range shown. The data shown in the background is the same as in Figure [5] The insets show the redshift distribution 
of LFs avalable in for the calculation in each band (i.e. Ao— AA < X obs < Ao+AA). 

despite these limitations, we obtain very comparable number counts to the ones in Section [4j agreeing to within 20% 
in the relevant magnitude range. 

B. CONSISTENCY NOTES 

K-correction 

Calculating the absolute magnitudes of a galaxy sample requires a k-correction to a c count for the offset in the 
rest-frame and the observed SED due to the cosmological redshift (e.g. iHogg et al.l (|2002f) ; iBlanton et al.l (|2003aD ) 

M x = m x - DM(z) - K{\ x ,\x>) (Bl) 

where X refers the band of interest. The k-correction can be written (in AB magnitudes) 

K{z) = {m x > - m x ) - 2.51og 10 (l + z) (B2) 

where mx is the observed brightness of a galaxy at redshift z and m' x is its rest-frame brightness in X-band. The exact 
value of the k-correction requires knowledge of the spectral energy distribution (SED) of the source and is commonly 
evaluated by assuming a template SEDs based on the galaxy type/color. This treatment is fairly reliable for low- 
z galaxies but the correction can become large for high-z galaxies and dominate the uncertainty in the derived LF, 
especially in the blue bands. Recent multiband photometric surveys offer a robust way of reducing this SED dependency 
by utilizing magnitudes in multiple bands to constrain the best-fit SED . Not only does multiband coverage indicate 
SED shape but when probing the LF in the rest- frame band Y centered at Xy-, the galaxy flux can be sampled in 
the band X which is closest to Ay(l + z). In other words, the observed filter (X) that best matches the redshifted 
rest-frame band of interest is the one that minimizes \Xx — Ay(l + z)\. The k-correction needed then becomes the 
matter of setting this quantity to exactly zero which is typically a small correction. We can rewrite Equation (|Blj) in 
this framework 

M Y = m x -DM(z)-K{\ x ,\y(1 + z)) (B3) 

where the SED dependence of the k-correction is now small even at high redshifts. Backtracking the original procedure 
to apparent magnitudes now requires simply K(z) — — 2.51og 10 (l + z) which we use in Equation[7] 

Photometric Systems 

Unfortunately, there is no photometric system which is universally accepted and the different ways used to evaluate 
the apparen t mag nitude of galaxies in the survey can introduce biases affecting the derived luminosity functions 
(see IBesselll (2005) for a review of photometric systems). As the flux from a galaxy diminishes from the center it 
will eventually drop below the background noise to be missed by the aperture. Photometric systems based on total 
magnitudes, such as Sersic, are usually preferred since they directly quantify the physical flux while apertures such as 
Kron and Petrosian will always suffer from missed light to some extent. However, t otal magnitudes t ypically assume 
an extrapolated profile which is mo del- dependent and has larger measurement errors ([Cole et al.ll200l . The Petrosian 
system can be advantageous since it compensat es for the effects of se eing by increasing the fraction of the light recovered 
from a galaxy when its angular siz e is small (IBlanton et al.l 120011). Despite this, Petrosian magnitudes are found to 
underestimate Sersic by 0.2 mag ([Strauss et al.l 120021 : IBlanton et aD 120011 ). Likewis e, 2MASS Krq n and isophotal 
magnitudes may a ccount for only 50-80% of the total flux in the most extreme cases (jAndreonl 120 02) ). For example, 
ISmith et all (|2009) show that their UKIDSS Petrosian magnitudes can be up to 0.5 mag fainter than 2MASS Kron 
magnitudes. The fraction of the lost flux increases towards faint er galaxies and m ay cause a systematic underestimation 
of the faint-end luminosities as well as the luminosity density. iHill et ail ([201 ll ) provide a good analysis of the effects 
of different photometric systems used in surveys. They find an overdensity of faint galaxies when compared with the 
best-fit Schechter function irrespective of the aperture system used and show that a Schechter function parametrization 



20 



does not provide a good fit at the faint-end. They also show that the use of a photometric systems based on total 
magnitudes (e.g. Sersic extrapolated) have a systematically steeper faint-end slope than photometric systems based 
on Kron or Petrosian magnitudes. They further show that the r-band Kron & Petrosian photometry u nderestimates 
the luminosity density by at least ~15% as they do not account for missing light. iBlanton et"al . (2003b!) show that the 
difference of the luminosity density resulting from Petrosian and Sersic magnitudes should be within <0.1 mag in the 
SDSS bands and not worth correcting for given the limitations of both systems. Still many authors apply a correction 
to estimate the total magnitudes i n ord e r to derive quan tities such as the luminosity density in physical units (e.g. 
iKochanek et~aTl (|200lD ; lBell et al.l (|2003ft ; lEke et all (|2005[ )V These can be as high as 0.3 mag in the K-band. It seems 
that uncertainties in the LF may be dominated by the aperture governing the fraction of flux recovered, especially at 
the faint end. 

Luminosity Function Estimators 

In this paper we use LFs derived from a variety of different LF estimators. The choice of LF estimator is unlikely 
to be a major source of discrepancy between the LFs derived by different authors although it can lead to different 
combinations of the Schechtcr parameters. The most commonly used met hods are i) the 1/Vm ax method (|Schmidtl 
1968), ii) the Sandage-Tammann-Yahil m aximum likelihood met hod (STY) flSandage et al.llT979h and hi) the Step Wise 
maximum Likelihood Method (SWLM) (lEfstathiou et al.lll988|) . The 1/V max method is reliable in the sense that it 
simultaneously gives the shape and normalization of the LF requiring no assumption on the parametric form for 
the LF. However, it suffers from systematic biases in the presence of density inhomogeneities in the observed field. 
The STY method is typically preferred when estimating the LF over mult iple fields since it has been shown to be 
unbiased to large scale structure and does not require binning of the data (Efstathi ou et al.l 119881 ). It does however 
require an assumption of a functional form of the luminosity function. The SWML metho d is widely use d since it 
makes no assumption of the LF shape while still being insensitive to large scale structure. IWillmerl (Jl997|) compare 
the properties of each LF estimator and show how different LF estimators tend to be biased towards the faint-end 
either overestimating or underestimating the slope, depending on the estimator and the unde rlying catalog. In orde r 
to minimize such e ffects one routinely co mpares the outputs of more than one method (e.g. iBouwens et al.l (2007); 
lllbert et al.l (120051 ): ICirasuolo et al.l (l20Tol) b 



